# Title: How and When Democratic Values Matter: Challenging the Effectiveness Centric Framework in Program Evaluation
# Yixin Liu, Heewon Lee, Frances Berry (LLB)
# Public Performance & Management Review
# Replication file

library(ggplot2)
library(lattice)
library(effsize)
library(effectsize)
library(stargazer)
library(Hmisc)
library(cregg)
library(ggpubr)
library(egg)
library(plyr)
library(dplyr)
library(multiwayvcov)
library(lmtest)
library(cowplot)
library(Rmisc)
library(yhat)

####Function####
element_textbox <- function(...) {
  el <- element_text(...)
  class(el) <- c("element_textbox", class(el))
  el
}

element_grob.element_textbox <- function(element, ...) {
  text_grob <- NextMethod()
  rect_grob <- element_grob(calc_element("strip.background", theme_bw()))
  
  ggplot2:::absoluteGrob(
    grid::gList(
      element_grob(calc_element("strip.background", theme_bw())),
      text_grob
    ),
    height = grid::grobHeight(text_grob), 
    width = grid::unit(1, "npc")
  )
}

# Setup working directory
WD = getwd()
setwd(WD)

####datefiles####
conjoint<-read.csv("Conjoint_LBB_PPMR.csv")
survey<-read.csv("Survey_LBB_PPMR.csv")

####clean data####
# drop manipulation check & attention test failure
conjoint<-subset(conjoint, mc.t==1 & at==1)
survey<-subset(survey, mc.t==1 & at==1)

# charaterize trust treatment variable
survey$rand<-ifelse(survey$rand==1,"LT Group","Control Group")

# generate binary variable for rating DV
conjoint$rdummy<-ifelse(conjoint$rate>50,1,0)

# convert to factor
conjoint$Inclusiveness[conjoint$Inclusiveness=="Diverse local communities"]<-"Local communities"
conjoint$Inclusiveness=as.factor(conjoint$Inclusiveness)
conjoint$Openness=as.factor(conjoint$Openness)
conjoint$Environmental.indicator=as.factor(conjoint$Environmental.indicator)
conjoint$Economic.indicator=as.factor(conjoint$Economic.indicator)

conjoint$prime1<-ifelse(conjoint$prime==1,"LT Group", "Control Group")
conjoint$prime1<-as.factor(conjoint$prime1)

# label attributes
label(conjoint$Inclusiveness)<-"Decision-making involves"
label(conjoint$Openness)<-"Implementation information is available to"
label(conjoint$Environmental.indicator)<-"Reduce annual CO2 emission"
label(conjoint$Economic.indicator)<-"Save schools’ annual expenses"

feature.labs <- c("Democratic Value: Openness", "Democratic Value: Inclusiveness", 
                  "Effectiveness: Environmental Performance", "Effectiveness: Economic Performance")
names(feature.labs) <- c("Implementation information is available to", 
                         "Decision-making involves",
                         "Reduce annual CO2 emission",
                         "Save schools’ annual expenses")

group_names <- c("Control Group","LT Group")

####Figure 2####

# Validity of the Low Trust Instrument
lt.df<-summarySE(survey, measurevar="trust", 
                 groupvars=c("rand"),
                 conf.interval = 0.95,
                 na.rm = T)

lt.plot<-ggplot(lt.df,aes(x=rand, y=trust, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(color=rand),size=5) +
  geom_pointrange(aes(color=rand, ymin=trust-ci, ymax=trust+ci)) +
  scale_color_manual(values = c("#00BFC4", "#F8766D"))+
  scale_y_continuous(lim=c(30,65)) +
  labs(x = element_blank(),
       y = "Degree of Trust",
       title = "Validity of the Low Trust Instrument")+
  annotate("label", x = 1.5, y = 62, label = "t-test, p-value = 0.00")+
  theme_article()
lt.plot<- lt.plot+theme(legend.position="none",
                        axis.text.x=element_text(size = 12),
                        plot.title = element_textbox(
                          hjust = 0.5, margin = margin(t = 5, b = 5)))
lt.plot

# Overall Probability of Program Support
trust.df<-summarySE(conjoint, measurevar="rdummy", 
                    groupvars=c("prime1"),
                    conf.interval = 0.95,
                    na.rm = T)

trust.rate.plot<-ggplot(trust.df,aes(x=prime1, y=rdummy, group=1)) +
  geom_line(linetype = "dashed") +
  geom_point(aes(color=prime1),size=5) +
  geom_pointrange(aes(color=prime1, ymin=rdummy-ci, ymax=rdummy+ci)) +
  scale_color_manual(values = c("#00BFC4", "#F8766D"))+
  labs(x = element_blank(),
       y = "Program Rating",
       title = "Overall Probability of Program Support")+
  annotate("label", x = 1.5, y = 0.82, label = "t-test, p-value = 0.00")+
  theme_article()
trust.rate.plot<- trust.rate.plot+theme(legend.position="none",
                                        axis.text.x=element_text(size = 12),
                                        plot.title = element_textbox(
                                          hjust = 0.5, margin = margin(t = 5, b = 5)))
trust.rate.plot

# combine
figure2<-ggarrange(lt.plot,trust.rate.plot,
                              ncol = 2,nrow = 1)
figure2<-plot_grid(figure2)
figure2
ggsave(file = "figure2.eps", width = 12, height = 4, dpi = 666)
ggsave(file = "figure2.jpg", width = 12, height = 4, dpi = 666)

####Figure 3####
# Change for Program Evaluation

# choice AMCE
choice.amce <- cj(conjoint, choice ~ Openness + Inclusiveness + 
                    Environmental.indicator + Economic.indicator, 
                  id = ~ IndID, estimate = "amce")

choice.plot<-plot(choice.amce, 
                  size = 4)+
  labs(title = "Probability of Program Support (Choice)",
       x = "Estimated AMCE")+
  theme(legend.position="none",
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  scale_color_manual(values = c("black","black","black","black"))+
  ggplot2::xlim(-0.3,0.3)

choice.plot<-choice.plot+theme(plot.title = element_textbox(
  hjust = 0.5, margin = margin(t = 5, b = 5)))

choice.plot

# rating AMCE
rate.amce <- cj(conjoint, rdummy ~ Openness + Inclusiveness + 
                  Environmental.indicator + Economic.indicator, 
                id = ~ IndID, estimate = "amce")

rate.plot<-plot(rate.amce, 
                size = 4)+
  labs(title = "Probability of Program Support (Rating)",
       x = "Estimated AMCE")+
  theme(legend.position="none",
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  scale_color_manual(values = c("black","black","black","black"))+
  ggplot2::xlim(-0.3,0.3)

rate.plot<-rate.plot+theme(plot.title = element_textbox(
  hjust = 0.5, margin = margin(t = 5, b = 5)))

rate.plot

# combine
figure3<-ggarrange(choice.plot,
                     rate.plot+theme(axis.text.y = element_blank(),
                                     axis.ticks.y = element_blank()),
                     ncol = 2,nrow = 1)
figure3<-plot_grid(figure3)
figure3
ggsave(file = "figure3.eps", width = 12, height = 7, dpi = 666)
ggsave(file = "figure31.jpg", width = 12, height = 7, dpi = 666)

####Figure 4####
# choice marginal means
choice.mm <- cj(conjoint, choice ~ Openness + Inclusiveness + 
                  Environmental.indicator + Economic.indicator, 
                id = ~ IndID, estimate = "mm", by = ~prime1)

choice.mm.plot<-plot(choice.mm, 
                     group = "prime1", 
                     size = 4)+
  labs(title = "Probability of Program Support (Choice)\nby Trust Levels",
       x = "Marginal Mean")+
  theme(legend.title=element_blank(),
        legend.position=c(0,1),
        legend.justification=c(0,1),
        legend.background = element_rect(size=0.5, linetype="solid", 
                                         colour ="black"),
        legend.key = element_blank(),
        legend.text = element_text(size = 11),
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11),
        plot.title = element_textbox(
          hjust = 0.5, margin = margin(t = 5, b = 5)))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  guides(color = guide_legend(reverse=T))+
  scale_color_manual("prime1",
                     breaks = group_names,
                     values = c("#00BFC4", "#F8766D"))

choice.mm.plot

# rating marginal means
rate.mm <- cj(conjoint, rdummy ~ Openness + Inclusiveness + 
                Environmental.indicator + Economic.indicator, 
              id = ~ IndID, estimate = "mm", by = ~prime1)

rate.mm.plot<-plot(rate.mm, 
                   group = "prime1", 
                   size = 4)+
  labs(title = "Probability of Program Support (Rating)\nby Trust Levels",
       x = "Marginal Mean")+
  theme(legend.title=element_blank(),
        legend.position=c(0,1),
        legend.justification=c(0,1),
        legend.background = element_rect(size=0.5, linetype="solid", 
                                         colour ="black"),
        legend.key = element_blank(),
        legend.text = element_text(size = 11),
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11),
        plot.title = element_textbox(
          hjust = 0.5, margin = margin(t = 5, b = 5)))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  guides(color = guide_legend(reverse=T))+
  scale_color_manual("prime1",
                     breaks = group_names,
                     values = c("#00BFC4", "#F8766D"))
rate.mm.plot

# choice diff-in-marginal means
choice.mmdiff <- cj(conjoint, choice ~ Openness + Inclusiveness + 
                      Environmental.indicator + Economic.indicator, 
                    id = ~ IndID, estimate = "mm_diff", by = ~prime1)

choice.mmdiff.plot<-plot(choice.mmdiff, 
                         group = "prime1", 
                         size = 4)+
  labs(title = "Difference in Probability of Program Support (Choice)\nby Trust Levels (Ref. = Control Group)",
       x = "Difference in Marginal Mean")+
  theme(legend.title=element_blank(),legend.position="none",
        legend.text = element_text(size = 11),
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11),
        plot.title = element_textbox(
          hjust = 0.5, margin = margin(t = 5, b = 5)))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  scale_color_manual("prime1",
                     breaks = group_names,
                     values = c("black", "black"))+
  ggplot2::xlim(-0.12,0.12)

choice.mmdiff.plot

# rating diff-in-marginal means
rate.mmdiff <- cj(conjoint, rdummy ~ Openness + Inclusiveness + 
                    Environmental.indicator + Economic.indicator, 
                  id = ~ IndID, estimate = "mm_diff", by = ~prime1)

rate.mmdiff.plot<-plot(rate.mmdiff, 
                       group = "prime1", 
                       size = 4)+
  labs(title ="Difference in Probability of Program Support (Rating)\nby Trust Levels (Ref. = Control Group)",
       x = "Difference in Marginal Mean")+
  theme(legend.title=element_blank(),legend.position="none",
        legend.text = element_text(size = 11),
        panel.grid.major = element_blank(),
        strip.text.x=element_text(size = 12),
        axis.text.y=element_text(size = 11),
        plot.title = element_textbox(
          hjust = 0.5, margin = margin(t = 5, b = 5)))+
  geom_errorbarh(aes_string(xmin = "lower", xmax = "upper"),
                 colour = "black",size=0.5,height = 0,
                 position = ggstance::position_dodgev(height = 0.75))+
  scale_color_manual("prime1",
                     breaks = group_names,
                     values = c("black", "black"))+
  ggplot2::xlim(-0.12,0.12)

rate.mmdiff.plot

# combine
figure4<-ggarrange(choice.mm.plot,
                   choice.mmdiff.plot+theme(axis.text.y = element_blank(),
                                            axis.ticks.y = element_blank()),
                   rate.mm.plot,
                   rate.mmdiff.plot+theme(axis.text.y = element_blank(),
                                          axis.ticks.y = element_blank()),
                   ncol = 2,nrow = 2)
figure4<-plot_grid(figure4)

figure4
ggsave(file = "figure4.eps", width = 12, height = 10, dpi = 666)
ggsave(file = "figure4.jpg", width = 13, height = 10, dpi = 666)

####Supplemental Information####
#####Appendix C####
vars<-data.frame(survey$female,
                 survey$white,survey$blk,survey$his,survey$asian,survey$other,
                 survey$ideo,
                 survey$age,
                 survey$education,
                 survey$income)

Mean=sapply(vars, mean, na.rm = TRUE)
SD=sapply(vars, sd, na.rm = TRUE)
Min=sapply(vars, min, na.rm = TRUE)
Max=sapply(vars, max, na.rm = TRUE)

# randomization check
female.aov <- aov(female ~ rand, data = survey)
summary(female.aov)
female<-summary(female.aov)[[1]][,5][1]

wht.aov <- aov(white ~ rand, data = survey)
summary(wht.aov)
wht<-summary(wht.aov)[[1]][,5][1]
blk.aov <- aov(blk ~ rand, data = survey)
summary(blk.aov)
blk<-summary(blk.aov)[[1]][,5][1]
his.aov <- aov(his ~ rand, data = survey)
summary(his.aov)
his<-summary(his.aov)[[1]][,5][1]
asian.aov <- aov(asian ~ rand, data = survey)
summary(asian.aov)
asian<-summary(asian.aov)[[1]][,5][1]
other.aov <- aov(other ~ rand, data = survey)
summary(other.aov)
other<-summary(other.aov)[[1]][,5][1]

ideo.aov <- aov(ideo ~ rand, data = survey)
summary(ideo.aov)
ideo<-summary(ideo.aov)[[1]][,5][1]

age.aov <- aov(age ~ rand, data = survey)
summary(age.aov)
age<-summary(age.aov)[[1]][,5][1]

educ.aov <- aov(education ~ rand, data = survey)
summary(educ.aov)
educ<-summary(educ.aov)[[1]][,5][1]

inc.aov <- aov(income ~ rand, data = survey)
summary(inc.aov)
inc<-summary(inc.aov)[[1]][,5][1]


# sample table
des_stat<-matrix(NA,length(Mean),5)
colnames(des_stat)<-c("Mean","SD","Min","Max","Randomization Check (P-value)")
rownames(des_stat)<-c("Female","White","Black","Hispanic","Asian","Other",
                      "Ideology","Age","Education","Income")
des_stat[,1]<-Mean
des_stat[,2]<-SD
des_stat[,3]<-Min
des_stat[,4]<-Max
des_stat[,5]<-c(female,wht,blk,his,asian,other,ideo,age,educ,inc)

# write table
write.table(des_stat, file = "appendixC", sep = ",", quote = FALSE, row.names = T)

#####AppendixD####
# Distribution of Trust in U.S. Local Governments
appendixD<-ggplot(survey,aes(x=trust,color=rand,fill=rand))+
  geom_histogram(position="dodge",bins=20)+
  labs(x="Degree of Trust",y="Count")+
  theme_article()+
  theme(legend.title=element_blank(),
        legend.text = element_text(size=10),
        legend.justification=c(0,1), 
        legend.position=c(0,1),
        legend.background=element_blank())+
  scale_color_manual(values=c("grey50", "black"))+
  scale_fill_manual(values=c("grey95", "grey70"))+
  guides(color = guide_legend(reverse=T),fill=guide_legend(reverse=T))
appendixD<-ggpar(appendixD,font.tickslab=c(10,"plain","black"))
appendixD
ggsave(file = "appendixD.eps", width = 6, height = 3, dpi = 666)
ggsave(file = "appendixD.jpg", width = 6, height = 3, dpi = 666)
